home *** CD-ROM | disk | FTP | other *** search
/ Pascal Super Library / Pascal Super Library (CW International)(1997).bin / MATH / COMPLX / COMPLEXO.PAS < prev    next >
Pascal/Delphi Source File  |  1992-01-15  |  23KB  |  729 lines
  1. {$N+,E+} {Use math coprocessor, if any, emulate otherwise.}
  2. UNIT ComplexOps;
  3.  
  4.  {This UNIT provides complex arithmetic and transcendental functions.
  5.  
  6.   (C) Copyright 1990, 1992, Earl F. Glynn, Overland Park, KS.  Compuserve 73257,3527.
  7.   All rights reserved.  This program may be freely distributed only for
  8.   non-commercial use.
  9.  
  10.   Some ideas in this UNIT were borrowed from "A Pascal Tool for Complex
  11.   Numbers", Journal of Pascal, Ada, & Modula-2, May/June 1985, pp. 23-29.
  12.   Many complex formulas were taken from Chapter 4, "Handbook of Mathematical
  13.   Functions" (Ninth Printing), Abramowitz and Stegun (editors), Dover, 1972.}
  14.  
  15. INTERFACE
  16.  
  17.   TYPE
  18.     RealType    = DOUBLE;
  19.     ComplexForm = (polar,rectangular);
  20.     Complex =
  21.       RECORD
  22.         CASE form:  ComplexForm OF
  23.           rectangular:  (x,y    :  RealType);  {z = x + y*i}
  24.           polar      :  (r,theta:  RealType);  {z = r*CIS(theta)}
  25.       END;                 {where CIS(theta) = COS(theta) + SIN(theta)*i}
  26.                            {      theta = -PI..PI (in canonical form)}
  27.  
  28.   CONST
  29.     MaxTerm    : BYTE     = 35;
  30.     EpsilonSqr : RealType = 1.0E-20;
  31.     Infinity   : RealType = 1.0E25;    {virtual infinity}
  32.  
  33.                              {complex number definition/conversion/output}
  34.   PROCEDURE CConvert (VAR z:  Complex; f:  ComplexForm);
  35.   PROCEDURE CSet (VAR z:  Complex; a,b:  RealType; f:  ComplexForm);
  36.   FUNCTION  CStr (z:  Complex; w,d:  BYTE; f:  ComplexForm):  STRING;
  37.  
  38.                              {complex arithmetic}
  39.   PROCEDURE CAdd  (VAR z:  Complex; a,b:  Complex);    {z = a + b}
  40.   PROCEDURE CDiv  (VAR z:  Complex; a,b:  Complex);    {z = a / b}
  41.   PROCEDURE CMult (VAR z:  Complex; a,b:  Complex);    {z = a * b}
  42.   PROCEDURE CSub  (VAR z:  Complex; a,b:  Complex);    {z = a - b}
  43.   PROCEDURE CNeg  (VAR z:  Complex; a  :  Complex);    {z = -a   }
  44.  
  45.                              {complex natural log, exponential}
  46.   PROCEDURE CLn   (VAR fn :  Complex; z:  Complex);   {fn  = ln(z) }
  47.   PROCEDURE CExp  (VAR z  :  Complex; a:  Complex);   {z   = exp(a)}
  48.   PROCEDURE CPwr  (VAR z  :  Complex; a,b:  Complex); {z   = a^b   }
  49.  
  50.                              {complex trig functions}
  51.   PROCEDURE CCos  (VAR z:  Complex; a:  Complex);     {z   = cos(a)}
  52.   PROCEDURE CSin  (VAR z:  Complex; a:  Complex);     {z   = sin(a)}
  53.   PROCEDURE CTan  (VAR z:  Complex; a:  Complex);     {z   = tan(a)}
  54.  
  55.   PROCEDURE CSec  (VAR z:  Complex; a:  Complex);     {z   = sec(a)}
  56.   PROCEDURE CCsc  (VAR z:  Complex; a:  Complex);     {z   = csc(a)}
  57.   PROCEDURE CCot  (VAR z:  Complex; a:  Complex);     {z   = cot(a)}
  58.  
  59.                              {complex hyperbolic functions}
  60.   PROCEDURE CCosh (VAR z:  Complex; a:  Complex);     {z   = cosh(a)}
  61.   PROCEDURE CSinh (VAR z:  Complex; a:  Complex);     {z   = sinh(a)}
  62.   PROCEDURE CTanh (VAR z:  Complex; a:  Complex);     {z   = tanh(a)}
  63.  
  64.   PROCEDURE CSech (VAR z:  Complex; a:  Complex);     {z   = sech(a)}
  65.   PROCEDURE CCsch (VAR z:  Complex; a:  Complex);     {z   = csch(a)}
  66.   PROCEDURE CCoth (VAR z:  Complex; a:  Complex);     {z   = coth(a)}
  67.  
  68.                              {miscellaneous complex functions}
  69.   FUNCTION  CAbs (z:  Complex):  RealType;                 {CAbs = |z|}
  70.   FUNCTION  CAbsSqr (z:  Complex):  RealType;           {CAbsSqr = |z|^2}
  71.   PROCEDURE CIntPwr (VAR z:  Complex; a:  Complex; n:  INTEGER); {z = a^n}
  72.   PROCEDURE CRealPwr (VAR z:  Complex; a:  Complex; x:  RealType); {z = a^x}
  73.   PROCEDURE CConjugate (VAR z:  Complex; a:  Complex);        {z = a*}
  74.   PROCEDURE CSqrt (VAR z:  Complex; a:  Complex);      {z = SQRT(a)}
  75.   PROCEDURE CRoot (VAR z:  Complex; a:  Complex; k,n:  WORD); {z = a^(1/n)}
  76.  
  77.                              {complex Bessel functions of order zero}
  78.   PROCEDURE CI0   (VAR sum:  Complex; z:  Complex);  {sum = I0(z)}
  79.   PROCEDURE CJ0   (VAR sum:  Complex; z:  Complex);  {sum = J0(z)}
  80.  
  81.   PROCEDURE CLnGamma (VAR z:  Complex; a:  Complex);
  82.   PROCEDURE CGamma   (VAR z:  Complex; a:  Complex);
  83.  
  84.                                   {treat "fuzz" of real numbers}
  85.   PROCEDURE CDeFuzz (VAR z:  Complex);
  86.   FUNCTION  DeFuzz (x:  RealType):  RealType;
  87.   PROCEDURE SetFuzz (value:  RealType);
  88.  
  89.                                   {miscellaneous}
  90.   FUNCTION FixAngle (theta:  RealType):  RealType;    {-PI < theta <= PI}
  91.   FUNCTION Pwr (x,y:  RealType):  RealType;    {Pwr = x^y}
  92.   FUNCTION Log10 (x:  RealType):  RealType;
  93.   FUNCTION LongMod (l1,l2:  LongInt):  LongInt;
  94.   FUNCTION Cosh (x:  RealType):  RealType;
  95.   FUNCTION Sinh (x:  RealType):  RealType;
  96.  
  97. IMPLEMENTATION
  98.  
  99.   VAR
  100.     fuzz     :  RealType;
  101.     Cone     :  Complex;
  102.     Cinfinity:  Complex;
  103.     Czero    :  Complex;
  104.     hLnPI    :  RealType;
  105.     hLn2PI   :  RealType;
  106.     ln2      :  RealType;
  107.  
  108.                              {complex number definition/conversion/output}
  109.   PROCEDURE CConvert (VAR z:  Complex; f:  ComplexForm);
  110.     VAR a:  Complex;
  111.   BEGIN
  112.     IF   z.form = f
  113.     THEN CDeFuzz (z)
  114.     ELSE BEGIN
  115.       CASE z.form OF
  116.         polar:            {polar-to-rectangular conversion}
  117.           BEGIN
  118.             a.form := rectangular;
  119.             a.x := z.r * COS(z.theta);
  120.             a.y := z.r * SIN(z.theta)
  121.           END;
  122.         rectangular:      {rectangular-to-polar conversion}
  123.           BEGIN
  124.             a.form := polar;
  125.             IF   DeFuzz(z.x) = 0.0
  126.             THEN BEGIN
  127.               IF   DeFuzz(z.y) = 0.0
  128.               THEN BEGIN
  129.                 a.r     := 0.0;
  130.                 a.theta := 0.0
  131.               END
  132.               ELSE
  133.                 IF   z.y > 0.0
  134.                 THEN BEGIN
  135.                   a.r := z.y;
  136.                   a.theta := 0.5*PI
  137.                 END
  138.                 ELSE BEGIN
  139.                   a.r := -z.y;
  140.                   a.theta := -0.5*PI
  141.                 END
  142.             END
  143.             ELSE BEGIN
  144.               a.r := CAbs(z);
  145.               a.theta := ARCTAN(z.y/z.x);   {4th/1st quadrant -PI/2..PI/2}
  146.               IF   z.x < 0.0                {2nd/3rd quadrants}
  147.               THEN
  148.                 IF   z.y >= 0.0
  149.                 THEN a.theta :=  PI + a.theta {2nd quadrant:  PI/2..PI}
  150.                 ELSE a.theta := -PI + a.theta {3rd quadrant: -PI..-PI/2}
  151.             END
  152.           END;
  153.       END;
  154.       CDeFuzz (a);
  155.       z := a
  156.     END
  157.   END {CConvert};
  158.  
  159.   PROCEDURE CSet (VAR z:  Complex; a,b:  RealType; f:  ComplexForm);
  160.   BEGIN
  161.     z.form := f;
  162.     CASE f OF
  163.       polar:
  164.         BEGIN
  165.           z.r := a;
  166.           z.theta := b
  167.         END;
  168.       rectangular:
  169.         BEGIN
  170.           z.x := a;
  171.           z.y := b
  172.         END;
  173.     END
  174.   END {CSet};
  175.  
  176.   FUNCTION  CStr (z:  Complex; w,d:  BYTE; f:  ComplexForm):  STRING;
  177.     VAR s1,s2:  STRING;
  178.   BEGIN
  179.     CConvert (z,f);
  180.     CASE f OF
  181.       polar:
  182.         BEGIN
  183.           Str (z.r:w:d, s1);
  184.           Str (z.theta:w:d, s2);
  185.           CStr := s1+'*CIS('+s2+')'
  186.         END;
  187.       rectangular:
  188.         BEGIN
  189.           Str (z.x:w:d, s1);
  190.           Str (ABS(z.y):w:d, s2);
  191.           IF   z.y >= 0
  192.           THEN CStr := s1+' +'+s2+'i'
  193.           ELSE CStr := s1+' -'+s2+'i'
  194.         END
  195.     END
  196.   END {CStr};
  197.  
  198.                                   {complex arithmetic}
  199.   PROCEDURE CAdd  (VAR z:  Complex; a,b:  Complex);    {z = a + b}
  200.   BEGIN                                               {complex addition}
  201.     CConvert (a,rectangular);
  202.     CConvert (b,rectangular);
  203.     z.form := rectangular;
  204.     z.x := a.x + b.x;   {real part}
  205.     z.y := a.y + b.y;   {imaginary part}
  206.   END {CAdd};
  207.  
  208.   PROCEDURE CDiv  (VAR z:  Complex; a,b:  Complex);    {z = a / b}
  209.     VAR temp:  RealType;
  210.   BEGIN
  211.     CConvert (b,a.form);    {arbitrarily convert one to type of other}
  212.     z.form := a.form;
  213.     CASE a.form OF
  214.       polar:
  215.         BEGIN
  216.           z.r := a.r / b.r;
  217.           z.theta := FixAngle(a.theta - b.theta)
  218.         END;
  219.       rectangular:
  220.         BEGIN
  221.           temp := SQR(b.x) + SQR(b.y);
  222.           z.x := (a.x*b.x + a.y*b.y) / temp;
  223.           z.y := (a.y*b.x - a.x*b.y) / temp
  224.         END
  225.     END
  226.   END {CDiv};
  227.  
  228.   PROCEDURE CMult (VAR z:  Complex; a,b:  Complex);    {z = a * b}
  229.   BEGIN
  230.     CConvert (b,a.form);    {arbitrarily convert one to type of other}
  231.     z.form := a.form;
  232.     CASE a.form OF
  233.       polar:
  234.         BEGIN
  235.           z.r := a.r * b.r;
  236.           z.theta := FixAngle(a.theta + b.theta)
  237.         END;
  238.       rectangular:
  239.         BEGIN
  240.           z.x := a.x*b.x - a.y*b.y;
  241.           z.y := a.x*b.y + a.y*b.x
  242.         END
  243.     END
  244.   END {CMult};
  245.  
  246.   PROCEDURE CSub  (VAR z:  Complex; a,b:  Complex);    {z = a - b}
  247.   BEGIN                                               {complex subtraction}
  248.     CConvert (a,rectangular);
  249.     CConvert (b,rectangular);
  250.     z.form := rectangular;
  251.     z.x := a.x - b.x;   {real part}
  252.     z.y := a.y - b.y;   {imaginary part}
  253.   END {CSub};
  254.  
  255.   PROCEDURE CNeg  (VAR z:  Complex; a  :  Complex);    {z = -a   }
  256.   BEGIN
  257.     z.form := a.form;
  258.     CASE a.form OF
  259.       polar:
  260.         BEGIN
  261.           z.r := a.r;
  262.           z.theta := FixAngle(a.theta + PI)
  263.         END;
  264.       rectangular:
  265.         BEGIN
  266.           z.x := -a.x;
  267.           z.y := -a.y
  268.         END
  269.     END
  270.   END {CNeg};
  271.                                   {complex natural log, exponential}
  272.   PROCEDURE CLn (VAR fn:  Complex; z:  Complex);  {fn  = ln(z)}
  273.   BEGIN  {Abramowitz formula 4.1.2 on p. 67}
  274.     CConvert (z,polar);
  275.     fn.form := rectangular;
  276.     fn.x := LN(z.r);
  277.     fn.y := FixAngle(z.theta)
  278.   END {CLn};  {principal value only}
  279.  
  280.   PROCEDURE CExp  (VAR z  :  Complex; a:  Complex);   {z   = exp(a)}
  281.     VAR
  282.       temp:  RealType;
  283.   BEGIN  {Euler's Formula:  Abramowitz formula 4.3.47 on p. 74}
  284.     CConvert (a,rectangular);
  285.     temp := EXP(a.x);
  286.     CSet (z, temp*COS(a.y),temp*SIN(a.y), rectangular)
  287.   END {CExp};
  288.  
  289.   PROCEDURE CPwr  (VAR z  :  Complex; a,b:  Complex); {z   = a^b   }
  290.     VAR
  291.       blna,lna:  Complex;
  292.   BEGIN  {Abramowitz formula 4.2.7 on p. 69}
  293.     CDeFuzz (a);
  294.     CDeFuzz (b);
  295.     IF   CAbsSqr(a) = 0.0
  296.     THEN
  297.       IF    (CAbsSqr(b) = 0.0)
  298.       THEN  z := Cone                   {lim a^a = 1 as a -> 0}
  299.       ELSE  z := Czero                  {0^b = 0, b <> 0}
  300.     ELSE BEGIN
  301.       CLn (lna,a);
  302.       CMult (blna,b,lna);
  303.       CExp (z, blna)
  304.     END
  305.   END {CPwr};
  306.                                   {complex trig functions}
  307.   PROCEDURE CCos  (VAR z:  Complex; a:  Complex);     {z   = cos(a)}
  308.   BEGIN  {Abramowitz formula 4.3.56 on p. 74}
  309.     CConvert (a,rectangular);
  310.     CSet (z, COS(a.x)*COSH(a.y), -SIN(a.x)*SINH(a.y), rectangular)
  311.   END {CCos};
  312.  
  313.   PROCEDURE CSin  (VAR z:  Complex; a:  Complex);     {z   = sin(a)}
  314.   BEGIN  {Abramowitz formula 4.3.55 on p. 74}
  315.     CConvert (a,rectangular);
  316.     CSet (z, SIN(a.x)*COSH(a.y), COS(a.x)*SINH(a.y), rectangular)
  317.   END {CSin};
  318.  
  319.   PROCEDURE CTan  (VAR z:  Complex; a:  Complex);     {z   = tan(a)}
  320.     VAR
  321.       temp:  RealType;
  322.   BEGIN  {Abramowitz formula 4.3.57 on p. 74}
  323.     CConvert (a,rectangular);
  324.     temp := COS(2.0*a.x) + COSH(2.0*a.y);
  325.     IF   DeFuzz(temp) <> 0.0
  326.     THEN BEGIN
  327.       CSet (z,SIN(2.0*a.x)/temp,SINH(2.0*a.y)/temp,rectangular)
  328.     END
  329.     ELSE z := Cinfinity
  330.   END {CTan};
  331.  
  332.   PROCEDURE CSec  (VAR z:  Complex; a:  Complex);     {z   = sec(a)}
  333.     VAR
  334.       temp:  Complex;
  335.   BEGIN  {Abramowitz formula 4.3.5 on p. 72}
  336.     CCos (temp, a);
  337.     IF   DeFuzz( Cabs(temp) ) <> 0.0
  338.     THEN CDiv (z, Cone,temp)
  339.     ELSE z := Cinfinity
  340.   END {CSec};
  341.  
  342.   PROCEDURE CCsc  (VAR z:  Complex; a:  Complex);     {z   = csc(a)}
  343.     VAR
  344.       temp:  Complex;
  345.   BEGIN  {Abramowitz formula 4.3.4 on p. 72}
  346.     CSin (temp, a);
  347.     IF   DeFuzz( Cabs(temp) ) <> 0.0
  348.     THEN CDiv (z, Cone,temp)
  349.     ELSE z := Cinfinity
  350.   END {CCsc};
  351.  
  352.   PROCEDURE CCot  (VAR z:  Complex; a:  Complex);     {z   = cot(a)}
  353.     VAR
  354.       temp:  RealType;
  355.   BEGIN  {Abramowitz formula 4.3.58 on p. 74}
  356.     CConvert (a,rectangular);
  357.     temp := COSH(2.0*a.y) - COS(2.0*a.x);
  358.     IF   DeFuzz(temp) <> 0.0
  359.     THEN BEGIN
  360.       CSet (z,SIN(2.0*a.x)/temp,-SINH(2.0*a.y)/temp,rectangular)
  361.     END
  362.     ELSE z := Cinfinity
  363.   END {CCot};
  364.  
  365.                                   {complex hyperbolic functions}
  366.   PROCEDURE CCosh (VAR z:  Complex; a:  Complex);     {z   = cosh(a)}
  367.   BEGIN  {Abramowitz formula 4.5.50 on p. 84}
  368.     CConvert (a,rectangular);
  369.     CSet (z, COSH(a.x)*COS(a.y), SINH(a.x)*SIN(a.y), rectangular)
  370.   END {CCosh};
  371.  
  372.   PROCEDURE CSinh (VAR z:  Complex; a:  Complex);     {z   = sinh(a)}
  373.   BEGIN  {Abramowitz formula 4.5.49 on p.84}
  374.     CConvert (a,rectangular);
  375.     CSet (z, SINH(a.x)*COS(a.y), COSH(a.x)*SIN(a.y), rectangular)
  376.   END {CSinh};
  377.  
  378.   PROCEDURE CTanh (VAR z:  Complex; a:  Complex);     {z   = tanh(a)}
  379.     VAR
  380.       temp:  RealType;
  381.   BEGIN  {Abramowitz formula 4.5.51 on p. 84}
  382.     CConvert (a,rectangular);
  383.     temp := COSH(2.0*a.x) + COS(2.0*a.y);
  384.     IF   DeFuzz(temp) <> 0.0
  385.     THEN BEGIN
  386.       CSet (z,SINH(2.0*a.x)/temp,SIN(2.0*a.y)/temp,rectangular)
  387.     END
  388.     ELSE z := Cinfinity
  389.   END {CTanh};
  390.  
  391.   PROCEDURE CSech (VAR z:  Complex; a:  Complex);     {z   = sech(a)}
  392.     VAR
  393.       temp:  Complex;
  394.   BEGIN  {Abramowitz formula 4.5.5 on p. 83}
  395.     CCosh (temp, a);
  396.     IF   DeFuzz( Cabs(temp) ) <> 0.0
  397.     THEN CDiv (z, Cone,temp)
  398.     ELSE z := Cinfinity
  399.   END {CSec};
  400.  
  401.   PROCEDURE CCsch (VAR z:  Complex; a:  Complex);     {z   = csch(a)}
  402.     VAR
  403.       temp:  Complex;
  404.   BEGIN  {Abramowitz formula 4.5.4 on p. 83}
  405.     CSinh (temp, a);
  406.     IF   DeFuzz( Cabs(temp) ) <> 0.0
  407.     THEN CDiv (z, Cone,temp)
  408.     ELSE z := Cinfinity
  409.   END {CCsch};
  410.  
  411.   PROCEDURE CCoth (VAR z:  Complex; a:  Complex);     {z   = coth(a)}
  412.     VAR
  413.       temp:  RealType;
  414.   BEGIN  {Abramowitz formula 4.5.52 on p. 84}
  415.     CConvert (a,rectangular);
  416.     temp := COSH(2.0*a.x) - COS(2.0*a.y);
  417.     IF   DeFuzz(temp) <> 0.0
  418.     THEN BEGIN
  419.       CSet (z,SINH(2.0*a.x)/temp,-SIN(2.0*a.y)/temp,rectangular)
  420.     END
  421.     ELSE z := Cinfinity
  422.   END {CCoth};
  423.  
  424.                                   {miscellaneous complex functions}
  425.   FUNCTION CAbs (z:  Complex):  RealType;                  {CAbs = |z|}
  426.   BEGIN
  427.     CASE z.form OF
  428.       rectangular:  CAbs := SQRT( SQR(z.x) + SQR(z.y) );
  429.       polar:        CAbs := z.r
  430.     END
  431.   END {CAbs};
  432.  
  433.   FUNCTION CAbsSqr (z:  Complex):  RealType;            {CAbsSqr = |z|^2}
  434.   BEGIN
  435.     CASE z.form OF
  436.       rectangular:  CAbsSqr := SQR(z.x) + SQR(z.y);
  437.       polar:        CAbsSqr := SQR(z.r)
  438.     END
  439.   END {CAbsSqr};
  440.  
  441.   PROCEDURE CIntPwr (VAR z:  Complex; a:  Complex; n:  INTEGER); {z = a^n}
  442.     {CIntPwr directly applies DeMoivre's theorem to calculate
  443.      an integer power of a complex number.  The formula holds
  444.      for both positive and negative values of 'n'.}
  445.   BEGIN
  446.     IF   CAbsSqr(a) = 0.0
  447.     THEN
  448.       IF    n = 0
  449.       THEN  z := Cone                   {lim a^a = 1 as a -> 0}
  450.       ELSE  z := Czero                  {0^n = 0, except for 0^0=1}
  451.     ELSE BEGIN
  452.       CConvert (a,polar);
  453.       z.form := polar;
  454.       z.r := Pwr(a.r,n);
  455.       z.theta := FixAngle(n*a.theta)
  456.     END
  457.   END {CIntPwr};
  458.  
  459.   PROCEDURE CRealPwr (VAR z:  Complex; a:  Complex; x:  RealType); {z = a^x}
  460.   BEGIN
  461.     IF   CAbsSqr(a) = 0.0
  462.     THEN
  463.       IF    Defuzz(x) = 0.0
  464.       THEN  z := Cone                   {lim a^a = 1 as a -> 0}
  465.       ELSE  z := Czero                  {0^n = 0, except for 0^0=1}
  466.     ELSE BEGIN
  467.       CConvert (a,polar);
  468.       z.form := polar;
  469.       z.r := Pwr(a.r,x);
  470.       z.theta := FixAngle(x*a.theta)
  471.     END
  472.   END {CRealPwr};
  473.  
  474.   PROCEDURE CConjugate (VAR z:  Complex; a:  Complex);      {z = a*}
  475.   BEGIN
  476.     z.form := a.form;
  477.     CASE a.form OF
  478.       polar:
  479.         BEGIN
  480.           z.r := a.r;
  481.           z.theta := FixAngle(-a.theta)
  482.         END;
  483.       rectangular:
  484.         BEGIN
  485.           z.x := a.x;
  486.           z.y := -a.y
  487.         END
  488.     END
  489.   END {CConjugate};
  490.  
  491.   PROCEDURE CSqrt (VAR z:  Complex; a:  Complex);      {z = SQRT(a)}
  492.   BEGIN
  493.     CRoot (z, a, 0,2)  {return only one of the two values}
  494.   END {CSqrt};
  495.                                            {z = a^(1/n), n > 0}
  496.   PROCEDURE CRoot (VAR z:  Complex; a:  Complex; k,n:  WORD);
  497.     {CRoot can calculate all 'n' roots of 'a' by varying 'k' from 0..n-1.}
  498.     {This is another application of DeMoivre's theorem.  See CIntPwr.}
  499.   BEGIN
  500.     IF   CAbs(a) = 0.0
  501.     THEN z := Czero                   {0^z = 0, except 0^0 is undefined}
  502.     ELSE BEGIN
  503.       CConvert (a,polar);
  504.       z.form := polar;
  505.       z.r := Pwr(a.r,1.0/n);
  506.       z.theta := FixAngle((a.theta + k*2.0*PI)/n)
  507.     END
  508.   END {CRoot};
  509.  
  510.                              {complex Bessel functions of order zero}
  511.   PROCEDURE CI0   (VAR sum:  Complex; z:  Complex);  {sum = I0(z)}
  512.     {I0(z) = Σ ( (¼z^2)^k / (k!)^2 ), k=0,1,2,...,∞}
  513.     VAR
  514.       i      :  BYTE;
  515.       SizeSqr:  RealType;
  516.       term   :  Complex;
  517.       zSQR25 :  Complex;
  518.   BEGIN
  519.     CConvert (z,rectangular);
  520.     sum := Cone;                       {term 0}
  521.     Cmult (zSQR25, z,z);
  522.     zSQR25.x := 0.25 * zSQR25.x;
  523.     zSQR25.y := 0.25 * zSQR25.y;       {¼z^2}
  524.     term := zSQR25;
  525.     CAdd (sum, sum,zSQR25);            {term 1}
  526.     i := 1;
  527.     REPEAT
  528.       CMult (term,zSQR25,term);
  529.       INC (i);
  530.       term.x := term.x / SQR(i);
  531.       term.y := term.y / SQR(i);
  532.       CAdd (sum, sum,term);       {sum := sum + term}
  533.       SizeSqr := SQR(term.x) + SQR(term.y)
  534.     UNTIL (i > MaxTerm) OR (SizeSqr < EpsilonSqr)
  535.   END {CI0};
  536.  
  537.   PROCEDURE CJ0   (VAR sum:  Complex; z:  Complex);  {sum = J0(z)}
  538.     {J0(z) = Σ ( (-1)^k * (¼z^2)^k / (k!)^2 ), k=0,1,2,...,∞}
  539.     VAR
  540.       addflag:  BOOLEAN;
  541.       i      :  BYTE;
  542.       SizeSqr:  RealType;
  543.       term   :  Complex;
  544.       zSQR25 :  Complex;
  545.   BEGIN
  546.     CConvert (z,rectangular);
  547.     sum := Cone;                       {term 0}
  548.     Cmult (zSQR25, z,z);
  549.     zSQR25.x := 0.25 * zSQR25.x;
  550.     zSQR25.y := 0.25 * zSQR25.y;       {¼z^2}
  551.     term := zSQR25;
  552.     CSub (sum, sum,zSQR25);            {term 1}
  553.     addflag := FALSE;
  554.     i := 1;
  555.     REPEAT
  556.       CMult (term,zSQR25,term);
  557.       INC (i);
  558.       addflag := NOT addflag;
  559.       term.x := term.x / SQR(i);
  560.       term.y := term.y / SQR(i);
  561.       IF   addflag
  562.       THEN CAdd (sum, sum,term)        {sum := sum + term}
  563.       ELSE CSub (sum, sum,term);       {sum := sum - term}
  564.       SizeSqr := SQR(term.x) + SQR(term.y)
  565.     UNTIL (i > MaxTerm) OR (SizeSqr < EpsilonSqr)
  566.   END {CJ0};
  567.  
  568.   PROCEDURE CApproxLnGamma (VAR sum:  Complex; z:  Complex);
  569.     {This is the approximation used in the National Bureau of
  570.      Standards "Table of the Gamma Function for Complex Arguments,"
  571.      Applied Mathematics Series 34, 1954.  The NBS table was created
  572.      using this approximation over the area  9 ≤ Re(z) ≤ 10 and
  573.      0 ≤ Im(z) ≤ 10.  Other table values were computed using the
  574.      relationship ln Γ(z+1) = ln z + ln Γ(z).}
  575.     CONST
  576.       c:  ARRAY[1..8] OF RealType
  577.        =  (1/12, -1/360, 1/1260, -1/1680, 1/1188, -691/360360,
  578.            1/156, -3617/122400);
  579.     VAR
  580.       i     :  WORD;
  581.       powers:  ARRAY[1..8] OF Complex;
  582.       temp1 :  Complex;
  583.       temp2 :  Complex;
  584.   BEGIN
  585.     CConvert (z,rectangular);
  586.     CLn  (temp1,z);                              {ln(z}
  587.     CSet (temp2, z.x-0.5, z.y, rectangular);     {z - 0.5}
  588.     CMult (sum, temp1,temp2);                    {(z - 0.5)*ln(z)}
  589.     CSub (sum, sum,z);                           {(z - 0.5)*ln(z) - z}
  590.     sum.x := sum.x + hLn2PI;
  591.  
  592.     temp1 := Cone;
  593.     CDiv (powers[1], temp1, z);                  {z^-1}
  594.     CMult (temp2, powers[1],powers[1]);          {z^-2}
  595.     FOR i := 2 TO 8 DO
  596.       CMult (powers[i], powers[i-1],temp2);
  597.     FOR i := 8 DOWNTO 1 DO BEGIN
  598.       CSet (temp1, c[i]*powers[i].x, c[i]*powers[i].y, rectangular);
  599.       CAdd (sum, sum,temp1);
  600.     END
  601.   END {CApproxLnGamma};
  602.  
  603.   PROCEDURE CLnGamma (VAR z:  Complex; a:  Complex);
  604.     VAR
  605.       lna :  Complex;
  606.       temp:  Complex;
  607.   BEGIN
  608.     CConvert (a, rectangular);
  609.  
  610.     IF   (a.x <= 0.0) AND (DeFuzz(a.y) = 0.0)
  611.     THEN
  612.       IF   DeFuzz(INT(a.x-1E-8) - a.x) = 0.0     {negative integer?}
  613.       THEN BEGIN
  614.         z := Cinfinity;
  615.         EXIT
  616.       END;
  617.  
  618.     IF   a.y < 0.0                     {3rd or 4th quadrant?}
  619.     THEN BEGIN
  620.       CConjugate (a, a);
  621.       CLnGamma (z, a);                 {try again in 1st or 2nd quadrant}
  622.       CConjugate (z, z)                {left this out!  1/3/91}
  623.     END
  624.     ELSE BEGIN
  625.       IF   a.x < 9.0                   {"left" of NBS table range}
  626.       THEN BEGIN
  627.         CLn (lna, a);
  628.         CSet (a, a.x+1.0, a.y, rectangular);
  629.         CLnGamma (temp,a);
  630.         CSub (z, temp,lna)
  631.       END
  632.       ELSE CApproxLnGamma (z,a)  {NBS table range:  9 ≤ Re(z) ≤ 10}
  633.     END
  634.   END {CLnGamma};
  635.  
  636.   PROCEDURE CGamma (VAR z:  Complex; a:  Complex);
  637.     VAR
  638.       lnz:  Complex;
  639.   BEGIN
  640.     CLnGamma (lnz,a);
  641.     IF   lnz.x >= 75.0       {arbitrary cutoff for infinity}
  642.     THEN z := Cinfinity
  643.     ELSE
  644.       IF   lnz.x < -200.0
  645.       THEN z := Czero        {avoid underflow}
  646.       ELSE CExp (z, lnz);
  647.   END {CGamma};
  648.  
  649.                                   {treat "fuzz" of real numbers}
  650.   PROCEDURE CDeFuzz (VAR z:  Complex);
  651.   BEGIN
  652.     CASE z.form OF
  653.       rectangular:
  654.         BEGIN
  655.           z.x := DeFuzz(z.x);
  656.           z.y := DeFuzz(z.y);
  657.         END;
  658.       polar:
  659.         BEGIN
  660.           z.r := DeFuzz(z.r);
  661.           IF   z.r = 0.0
  662.           THEN z.theta := 0.0     {canonical form when radius is zero}
  663.           ELSE z.theta := DeFuzz(z.theta)
  664.         END
  665.     END
  666.   END {CDeFuzz};
  667.  
  668.   FUNCTION  DeFuzz (x:  RealType):  RealType;
  669.   BEGIN
  670.     IF   ABS(x) < fuzz
  671.     THEN Defuzz := 0
  672.     ELSE Defuzz := x
  673.   END {Defuzz};
  674.  
  675.   PROCEDURE SetFuzz (value:  RealType);
  676.   BEGIN
  677.     fuzz := value
  678.   END {SetFuzz};
  679.  
  680.                                   {miscellaneous}
  681.   FUNCTION FixAngle (theta:  RealType):  RealType;    {-PI < theta <= PI}
  682.   BEGIN
  683.     WHILE theta > PI DO
  684.       theta := theta - 2.0*PI;
  685.     WHILE theta <= -PI DO
  686.       theta := theta + 2.0*PI;
  687.     FixAngle := DeFuzz(theta)
  688.   END {FixAngle};
  689.  
  690.   FUNCTION Pwr (x,y:  RealType):  RealType;        {Pwr = x^y}
  691.   BEGIN
  692.     IF   DeFuzz(x) = 0.0
  693.     THEN
  694.       IF   DeFuzz(y) = 0.0
  695.       THEN Pwr := 1.0    {0^0 = 1 (i.e., lim x^x = 1 as x -> 0)}
  696.       ELSE Pwr := 0.0
  697.     ELSE Pwr := EXP( LN(x)*y )
  698.   END {Pwr};
  699.  
  700.   FUNCTION Log10 (x:  RealType):  RealType;
  701.   BEGIN
  702.     Log10 := LN(x) / LN(10)
  703.   END {Log10};
  704.  
  705.   FUNCTION LongMod (l1,l2:  LongInt):  LongInt;
  706.   BEGIN
  707.     LongMod := l1 - l2*(l1 DIV l2)
  708.   END {LongMod};
  709.  
  710.   FUNCTION Cosh (x:  RealType):  RealType;
  711.   BEGIN
  712.     Cosh := 0.5*( EXP(x) + EXP(-x) )
  713.   END {Cosh};
  714.  
  715.   FUNCTION Sinh (x:  RealType):  RealType;
  716.   BEGIN
  717.     Sinh := 0.5*( EXP(x) - EXP(-x) )
  718.   END {Sinh};
  719.  
  720. BEGIN
  721.   SetFuzz (1.0E-12);
  722.   CSet ( Cone, 1.0, 0.0, rectangular);
  723.   CSet (Czero, 0.0, 0.0, rectangular);
  724.   CSet (Cinfinity, Infinity, Infinity, rectangular);
  725.   hLnPI := 0.5*(LN(PI));
  726.   hLn2PI := 0.5*(LN(2.0*PI));
  727.   ln2 := LN(2.0)
  728. END.
  729.